#include <feel/options.hpp>
#include <feel/feelalg/backend.hpp>
#include <feel/feeldiscr/functionspace.hpp>
#include <feel/feeldiscr/region.hpp>
#include <feel/feelpoly/im.hpp>
#include <feel/feelfilters/gmsh.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelpoly/polynomialset.hpp>
#include <feel/feelvf/vf.hpp>

#include <feel/feelfilters/geotool.hpp>
#include <feel/feeldiscr/operatorinterpolation.hpp>
/** use Feel namespace */
using namespace Feel;
using namespace Feel::vf;


namespace schwarz3d_cas2
{

inline
po::options_description
makeOptions()
{
  po::options_description code1doptions("Code1D options");
  code1doptions.add_options()
    ("hsize", po::value<double>()->default_value( 0.1 ), "mesh size")
    ("hsize1", po::value<double>()->default_value( 0.5 ), "mesh size1")
    ("hsize2", po::value<double>()->default_value( 0.4 ), "mesh size2")
    ("nu", po::value<double>()->default_value( 0.5 ), "nu")
    ("alphaX", po::value<double>()->default_value( 0.01 ), "alphaX")
    ("alphaY", po::value<double>()->default_value( -2.0 ), "alphaY")
    ("bccoeff", Feel::po::value<double>()->default_value( 100.0 ), "coeff for weak Dirichlet conditions")
    ("toler", Feel::po::value<double>()->default_value( 1e-6 ), "error tol")
    ("xc", Feel::po::value<double>()->default_value( 0.0 ), "param de recouvrement des maillages")
    ("penaldir", Feel::po::value<double>()->default_value( 10 ),
     "penalisation parameter for the weak boundary Dirichlet formulation")
    ("gmsh", po::value<float>()->default_value( 2.1 ), " version of gmsh( 2.0 or 2.1)")
    ;
  return code1doptions.add( Feel::feel_options() );
}

/**
 * This routine defines some information about the application like
 * authors, version, or name of the application. The data returned is
 * typically used as an argument of a Feel::Application subclass.
 *
 * \return some data about the application.
 */
inline
AboutData
makeAbout()
{
  AboutData about( "SchwarzMethod2" ,
		   "SchwarzMethod2" ,
		   "0.2",
		   "Code SchwarzMethod",
		   Feel::AboutData::License_GPL,
		   "Copyright (c) 2010 Universite Joseph Fourier");

  about.addAuthor("Vincent Chabannes", "developer", "vincent.chabannes@imag.fr", "");
  return about;

}


template<uint Dim,uint N,uint OrderGeo>
class Localisation
  :
  public Application
{
  typedef Application super;
public:

  //! Polynomial order \f$P_2\f$
  //    static const uint16_type Order = 5;
  ///static const uint Order = N;

  typedef double value_type;

  //---------------------------------------------------------------------------------------//

  typedef Backend<value_type> backend_type;
  typedef boost::shared_ptr<backend_type> backend_ptrtype;

  typedef backend_type::sparse_matrix_type sparse_matrix_type;
  typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
  typedef backend_type::vector_type vector_type;
  typedef backend_type::vector_ptrtype vector_ptrtype;

  //---------------------------------------------------------------------------------------//

  typedef Simplex<Dim,1,Dim> convex_1_type;
  typedef Mesh<convex_1_type> mesh_1_type;
  typedef boost::shared_ptr<mesh_1_type> mesh_1_ptrtype;

  typedef Simplex<Dim,OrderGeo,Dim> convex_N_type;
  typedef Mesh<convex_N_type> mesh_N_type;
  typedef boost::shared_ptr<mesh_N_type> mesh_N_ptrtype;

  //---------------------------------------------------------------------------------------//

    typedef Lagrange<N,Scalar,Continuous,PointSetFekete> basis_u_type;
  typedef bases<basis_u_type> basis_type;

  typedef FunctionSpace<mesh_1_type, basis_type> space_1_type;
  typedef FunctionSpace<mesh_N_type, basis_type> space_N_type;

  typedef boost::shared_ptr<space_1_type> space_1_ptrtype;
  typedef boost::shared_ptr<space_N_type> space_N_ptrtype;

  typedef typename space_1_type::element_type element_1_type;
  typedef typename space_N_type::element_type element_N_type;

  //---------------------------------------------------------------------------------------//

  typedef Exporter<mesh_1_type> export_1_type;
  typedef boost::shared_ptr<export_1_type> export_1_ptrtype;
  typedef Exporter<mesh_N_type> export_N_type;
  typedef boost::shared_ptr<export_N_type> export_N_ptrtype;

  /**
   * Constructor
   */
  Localisation( int argc, char** argv, AboutData const& ad,
		po::options_description const& od )
    :
    super( argc, argv, ad, od ),
    M_backend( backend_type::build( this->vm() ) ),
    M_backend2( backend_type::build( this->vm() ) ),
    meshSize( this->vm()["hsize"].template as<double>() ),
    meshSize1( this->vm()["hsize1"].template as<double>() ),
    meshSize2( this->vm()["hsize2"].template as<double>() ),
    exporter1( Exporter<mesh_N_type>::New( this->vm(), "domain_1" ) ),
    exporter2( Exporter<mesh_N_type>::New( this->vm(), "domain_2" ) )
  {
  }

  mesh_N_ptrtype createMesh( std::string __name,std::string __str) ;

  mesh_N_ptrtype createMesh( double meshSize );
  mesh_1_ptrtype createMeshO1( std::string __name,std::string __str );

  mesh_N_ptrtype createMesh1( double meshSize );
  mesh_N_ptrtype createMesh2( double meshSize );

    
  // run the application
  void runMultiDomain();

private:

  template<typename element_type, typename RhsExprType>
  void localProblem(element_type& U1 , RhsExprType );

  //void problem2(element_1_type& U2 , element_N_type& U1, element_N_type& lambda );

  void exportResults( element_N_type& U1,element_1_type& U2,double time );

private:

  // linear algebra backend
  backend_ptrtype M_backend,M_backend2;
  sparse_matrix_ptrtype A1;
  vector_ptrtype B1;
  sparse_matrix_ptrtype A2;
  vector_ptrtype B2;
  // mesh characteristic size
  double meshSize,meshSize1,meshSize2;

  export_N_ptrtype exporter1;
  export_N_ptrtype exporter2;
}; // Interpolate



template<uint Dim,uint N,uint OrderGeo>
typename Localisation<Dim,N,OrderGeo>::mesh_N_ptrtype
Localisation<Dim,N,OrderGeo>::createMesh( std::string __name,std::string __str)
{
  //typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
  //typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
  mesh_N_ptrtype mesh( new mesh_N_type );

  Gmsh gmsh;
  gmsh.setOrder(OrderGeo);

  std::string fname = gmsh.generate( __name, __str,false,false,false  );

  ImporterGmsh<mesh_N_type> import( fname );

  float gmsh_version = 2.1;//this->vm()["gmsh"].template as<float>();
  if (gmsh_version==2.0) import.setVersion( "2.0" );
  else if (gmsh_version==2.1) import.setVersion( "2.1" );

  mesh->accept( import );
  mesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
  mesh->updateForUse();

  return mesh;

}

template<uint Dim,uint N,uint OrderGeo>
typename Localisation<Dim,N,OrderGeo>::mesh_1_ptrtype
Localisation<Dim,N,OrderGeo>::createMeshO1( std::string __name,std::string __str)
{
  //typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
  //typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
  mesh_1_ptrtype mesh( new mesh_1_type );

  Gmsh gmsh;
  gmsh.setOrder(1);

  std::string fname = gmsh.generate( __name, __str ,false,false,false );

  ImporterGmsh<mesh_1_type> import( fname );

  float gmsh_version = 2.1;//this->vm()["gmsh"].template as<float>();
  if (gmsh_version==2.0) import.setVersion( "2.0" );
  else if (gmsh_version==2.1) import.setVersion( "2.1" );

  mesh->accept( import );
  mesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
  mesh->updateForUse();

  return mesh;

}


template<uint Dim,uint Order,uint OrderGeo>
typename Localisation<Dim,Order,OrderGeo>::mesh_N_ptrtype
Localisation<Dim,Order,OrderGeo>::createMesh1( double meshSize )
{
#if 0
    //double z=0.3;

        GeoTool::Node x1(-1,-1,-1);
        GeoTool::Node x2( 1,-1,-1);
        GeoTool::Node x3( 1, 1,-1);
        GeoTool::Node x4(-1, 1,-1);
        GeoTool::Node x5(-1,-1,1);
        GeoTool::Node x6( 1,-1,1);
        GeoTool::Node x7( 1, 1,1);
        GeoTool::Node x8(-1, 1,1);
        GeoTool::Hexaedre H(meshSize1,"Hexaedre",x1,x2,x3,x4,x5,x6,x7,x8);
        H.setMarker(_type="line",_name="Arete",_markerAll=true);
        H.setMarker(_type="surface",_name="Interface",_marker3=true);
        H.setMarker(_type="surface",_name="Bord",
                    _marker4=true,
                    _marker1=true,
                    _marker2=true,
                    _marker5=true,
                    _marker6=true);
        H.setMarker(_type="volume",_name="Omega1",_marker1=true);

        auto mesh = H.createMesh<mesh_N_type>("Omega1" + mesh_N_type::shape_type::name() );
        /*
        mesh_N_ptrtype mesh( createMesh( "Omega1" + mesh_N_type::shape_type::name(),
                                                     H.geoStr()
                                                      ));
        */
#else
  std::ostringstream ostr;

  ostr << "Mesh.MshFileVersion = 2;\n"
       << "Mesh.CharacteristicLengthExtendFromBoundary=1;\n"
       << "Mesh.CharacteristicLengthFromPoints=1;\n"
       << "Mesh.ElementOrder=1;\n"
       << "Mesh.SecondOrderIncomplete = 0;\n"
       << "Mesh.Algorithm = 6;\n"
       << "httt=" << meshSize1 <<";\n"       //h=0.5;
       << "l=0.75;\n"
       << "Point(1) = {-1,-1,-2, httt};\n"
       << "Point(2) = {0,-1,-2, httt};\n"
       << "Point(3) = {0,1,-2, httt};\n"
       << "Point(4) = {-1,1,-2, httt};\n"
       << "Point(5) = {-1,-1,2, httt};\n"
       << "Point(6) = {0,-1,2, httt};\n"
       << "Point(7) = {0,1,2, httt};\n"
       << "Point(8) = {-1,1,2, httt};\n"
       << "Point(9) = {-1,0,0, httt};\n"
       << "Point(10) = {-1,l,0, httt};\n"
       << "Point(11) = {-1,-l,0, httt};\n"
       << "Point(12) = {-1,0,l, httt};\n"
       << "Point(13) = {-1,0,-l, httt};\n"
       << "Line(1) = {1,2};\n"
       << "Line(2) = {2,3};\n"
       << "Line(3) = {3,4};\n"
       << "Line(4) = {4,1};\n"
       << "Line(5) = {5,6};\n"
       << "Line(6) = {6,7};\n"
       << "Line(7) = {7,8};\n"
       << "Line(8) = {8,5};\n"
       << "Line(9) = {1,5};\n"
       << "Line(10) = {2,6};\n"
       << "Line(11) = {3,7};\n"
       << "Line(12) = {4,8};\n"
       << "Circle(13) = {10, 9, 12};\n"
       << "Circle(14) = {12, 9, 11};\n"
       << "Circle(15) = {11, 9, 13};\n"
       << "Circle(16) = {13, 9, 10};\n"
       << "Line Loop(17) = {14, 15, 16, 13};\n"
       << "Plane Surface(18) = {17};\n"
       << "Line Loop(19) = {8, -9, -4, 12};\n"
       << "Plane Surface(20) = {17, 19};\n"
       << "Line Loop(21) = {7, -12, -3, 11};\n"
       << "Plane Surface(22) = {21};\n"
       << "Line Loop(23) = {11, -6, -10, 2};\n"
       << "Plane Surface(24) = {23};\n"
       << "Line Loop(25) = {5, -10, -1, 9};\n"
       << "Plane Surface(26) = {25};\n"
       << "Line Loop(27) = {5, 6, 7, 8};\n"
       << "Plane Surface(28) = {27};\n"
       << "Line Loop(29) = {1, 2, 3, 4};\n"
       << "Plane Surface(30) = {29};\n"
       << "Surface Loop(31) = {18, 20, 22, 28, 26, 24, 30};\n"
       << "Volume(32) = {31};\n"
       << "Physical Surface(\"Interface\") = {18};\n"
      //<< "Physical Surface(\"Entree\") = {28, 30};\n"
       << "Physical Surface(\"Bord\") = {  28, 30,22, 24, 26, 20};\n"
       << "Physical Volume(10) = {32};\n";
  std::cout << ostr.str();
  mesh_N_ptrtype mesh( createMesh( "Omega1" + mesh_N_type::shape_type::name(),ostr.str() ) );

#endif
  return mesh;


}
//----------------------------------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------------------------------//
//----------------------------------------------------------------------------------------------------------//

template<uint Dim,uint Order,uint OrderGeo>
typename Localisation<Dim,Order,OrderGeo>::mesh_N_ptrtype
Localisation<Dim,Order,OrderGeo>::createMesh2( double meshSize )
{
#if 0
    std::cout <<"\nOOLLAAA\n";
    //double z=0.3;

        GeoTool::Node centre(0,1,0);
        GeoTool::Node dir( 0,1,0);
        GeoTool::Node r(1.0);
        GeoTool::Node lg(4.0);
        GeoTool::Cylindre C(meshSize2,"Cylindre",centre,dir,r,lg);
        //C.setMarker(_type="line",_name="Arete",_markerAll=true);
        C.setMarker(_type="surface",_name="Interface",_marker1=true);
        C.setMarker(_type="surface",_name="Bord",
                    _marker2=true,
                    _marker3=true);
        C.setMarker(_type="volume",_name="Omega2",_marker1=true);
    std::cout <<"\nOOLLAAA01\n";

        auto mesh = C.createMesh<mesh_N_type>("Omega2" + mesh_N_type::shape_type::name() );
    std::cout <<"\nOOLLAAA2\n";

        /*
        mesh_N_ptrtype mesh( createMesh( "Omega2" + mesh_N_type::shape_type::name(),
                                                     H.geoStr()
                                                      ));
        */
#else
  std::ostringstream ostr;

  ostr << "Mesh.MshFileVersion = 2;\n"
       << "Mesh.CharacteristicLengthExtendFromBoundary=1;\n"
       << "Mesh.CharacteristicLengthFromPoints=1;\n"
       << "Mesh.ElementOrder=1;\n"
       << "Mesh.SecondOrderIncomplete = 0;\n"
       << "Mesh.Algorithm = 6;\n"
       << "lc = "<<meshSize2/*0.1*/<<";\n"
       << "lx = 0.75;\n"
       << "lz = 1.0;\n"
       << "delta_gap = 4.0;\n"
       << "Point(14) = {-1, 0, 0, lc};\n"
       << "Point(15) = {-1,-lx, 0, lc};\n"
       << "Point(16) = {-1, 0, -lx, lc};\n"
       << "Point(17) = {-1,lx, 0, lc};\n"
       << "Point(18) = {-1, 0, lx, lc};\n"
       << "Circle(17) = {15,14,16};\n"
       << "Circle(18) = {16,14,17};\n"
       << "Circle(19) = {17,14,18};\n"
       << "Circle(20) = {18,14,15};\n"
       << "Line Loop(32) = {19,20,17,18};\n"
       << "Plane Surface(60) = {32};\n"
       << "out23[] = Extrude {-delta_gap,0,0} { Surface{60}; } ;\n"
       << "Surface Loop(83) = {69, 60, 73, 77, 81, 82};\n"
       << "Volume(84) = {83};\n"
       << "Physical Surface(\"Interface\") = {60};\n"
      //<< "Physical Surface(\"Sortie\") = {out23[0]};\n"
       << "Physical Surface(\"Bord\") = {out23[0],out23[2],out23[3],out23[4],out23[5]};\n"
       << "Physical Volume(10) = {84};\n";

  mesh_N_ptrtype mesh( createMesh( "Omega2" + mesh_N_type::shape_type::name(),ostr.str() ) );

#endif
        return mesh;
}

template<uint Dim,uint Order,uint OrderGeo>
void
Localisation<Dim,Order,OrderGeo>::exportResults( element_N_type& U1, element_1_type& U2, double time  )
{
    //if ( exporter1->doExport() )
    {
      exporter1->step( time )->setMesh( U1.functionSpace()->mesh() );
      exporter1->step( time )->add( "u1", U1 );
      exporter1->save();

      exporter2->step( time )->setMesh( U2.functionSpace()->mesh() );
      exporter2->step( time )->add( "u2", U2 );
      exporter2->save();
    }
} // Test_Schwarz::exportResults




template<uint Dim,uint Order,uint OrderGeo>
void
Localisation<Dim,Order,OrderGeo>::runMultiDomain()
{

  if ( this->vm().count( "help" ) )
    {
      std::cout << this->optionsDescription() << "\n";
      return;
    }

  using namespace Feel::vf;

  //-----------------------------------------------------------------------------------//

  this->changeRepository( boost::format( "%1%/%2%/P%3%/h_%4%/" )
			  % this->about().appName()
			  % convex_N_type::name()
			  % Order
			  % this->vm()["hsize"].template as<double>()
			  );

  //-----------------------------------------------------------------------------------//

  mesh_N_ptrtype mesh1 = createMesh1( meshSize1 );
  mesh_N_ptrtype mesh2 = createMesh2( meshSize2 );

  space_N_ptrtype Xh1 = space_N_type::New( mesh1 );
  space_N_ptrtype Xh2 = space_N_type::New( mesh2 );

  element_N_type u1( Xh1, "U1" );
  element_N_type u2( Xh2, "U2" );
  element_N_type lambda( Xh1, "lambda" );
  element_N_type u2bis( Xh1, "u2bis" );

  element_N_type u1old( Xh1, "u1old" );
  element_N_type u2old( Xh2, "u2old" );


  //-----------------------------------------------------------------------------------//

  auto SolAnalytique = cos(M_PI*Px()/2.0)*cos(M_PI*Py()/2.0);


  typedef Backend<double> backend_type;
  typedef boost::shared_ptr<backend_type> backend_ptrtype;

  Feel::po::options_description meshHighOrderoptions("TestALE options");
  meshHighOrderoptions.add( Feel::feel_options() );
  po::variables_map vm;
  int ac = 0;
  char **av;
  po::store(po::parse_command_line(ac, av, meshHighOrderoptions), vm);

  backend_ptrtype backend( Backend<double>::build( /*BACKEND_PETSC*/ vm ) );

  //OperatorInterpolation<space_1_type, space_N_type> I( Xh2,Xh1 , backend );
  OperatorInterpolation<space_N_type, space_N_type> I( Xh2,Xh1 , backend );

  double erreur=1,errold=2;
  double err_tol=this->vm()["toler"].template as<double>();
  double err1,err2,norm1,norm2;
  int cpt=0;

  while (erreur>=err_tol && cpt <50)
    {
        errold=erreur;
      u1old=u1;
      u2old=u2;

      this->localProblem(u1,-idv(lambda) );

      //norm1 =math::sqrt( integrate( elements(mesh1), _Q<2*(Order+OrderGeo-1)>(), idv(u1)*idv(u1) ).evaluate()(0,0) );
      //err1=math::sqrt( integrate(elements(mesh1), _Q<2*(Order+OrderGeo-1)>(), (idv(u1)-idv(u1old))*(idv(u1)-idv(u1old)) ).evaluate()(0,0) );
      norm1 =math::sqrt( integrate( markedfaces(mesh1,mesh1->markerName("Interface")), _Q<2*(Order+OrderGeo-1)>(), idv(u1)*idv(u1) ).evaluate()(0,0) );
      err1=math::sqrt( integrate(markedfaces(mesh1,mesh1->markerName("Interface")), _Q<2*(Order+OrderGeo-1)>(), (idv(u1)-idv(u1old))*(idv(u1)-idv(u1old)) ).evaluate()(0,0) );

      if (norm1!=0)  err1=err1/norm1;

      this->localProblem(u2,idv(lambda));

      //norm2 =math::sqrt( integrate( elements(mesh2), _Q<2*Order>(), idv(u2)*idv(u2) ).evaluate()(0,0) );
      //err2=math::sqrt( integrate(elements(mesh2), _Q<2*Order>(), (idv(u2)-idv(u2old))*(idv(u2)-idv(u2old)) ).evaluate()(0,0) );
      norm2 =math::sqrt( integrate( markedfaces(mesh2,mesh2->markerName("Interface")), _Q<2*Order>(), idv(u2)*idv(u2) ).evaluate()(0,0) );
      err2=math::sqrt( integrate( markedfaces(mesh2,mesh2->markerName("Interface")), _Q<2*Order>(), (idv(u2)-idv(u2old))*(idv(u2)-idv(u2old)) ).evaluate()(0,0) );
      if (norm2!=0)  err2=err2/norm2;

      erreur=err1+err2;
      std::cout<<"\n erreur = "<<erreur<<" ; err1 = "<<err1<<" ; err2 = "<<err2;

      exportResults( u1, u2, cpt);
      if (erreur>errold ) erreur=1e-6;
      I.apply(u2,u2bis);

      //form1(Xh1,lambda.container())+= (0.5*idv(u1));

      lambda.add(0.5, u1);
      lambda.add(-0.5, u2bis);


      ++cpt;

    }

  //exportResults( u1, u2, 0);

}


//--//////--//////--//////--//////-------------------------------------------------------------------------------------------------//
//-//--//--//--//--//--//--//--//---------------------------------------------------------------------------------------//
////--//////--//////--//////--////------------------------------------------------------------------------------------------------//
template<uint Dim,uint Order,uint OrderGeo>
template<typename element_type,typename RhsExprType>
void
Localisation<Dim,Order,OrderGeo>::localProblem(element_type& u1,RhsExprType RhsExp ) {

  using namespace Feel::vf;

  //----------------------------------------------------------------------------------------------------------//

  typedef typename element_type::functionspace_type functionspace_type;
  typedef typename element_type::functionspace_ptrtype functionspace_ptrtype;
  typedef typename functionspace_type::mesh_type mesh_type;
  typedef typename functionspace_type::mesh_ptrtype mesh_ptrtype;

  functionspace_ptrtype Xh=u1.functionSpace();
  mesh_ptrtype mesh=Xh->mesh();

  const uint16_type GeoOrder = mesh_type::nOrder;

  element_type v(Xh,"v");


  auto pi = M_PI;
  //auto f = sin( pi*(Px()+1)/2.)*sin( pi*(Py()+1)/2.)*sin( pi*(Pz()+1)/2.);
  auto f= cst(1.);

  //----------------------------------------------------------------------------------------------------------//

  sparse_matrix_ptrtype A( M_backend->newMatrix(Xh,Xh) );

  form2( Xh, Xh, A, _init=true ) =
    integrate( elements(mesh), _Q<2*(Order+GeoOrder-1)>(),
	       gradt(u1)*trans(grad(v)) );

    A->close();

 //----------------------------------------------------------------------------------------------------------//

  vector_ptrtype B( M_backend->newVector( Xh ) );

  form1( Xh, B, _init=true ) =
    integrate( elements(mesh), f*id(v) );

  form1( Xh, B ) +=
    integrate( markedfaces(mesh,mesh->markerName("Interface")),
		 RhsExp*id(v) ); //-idv(lambda)*id(v) );

  B->close();

 //----------------------------------------------------------------------------------------------------------//

    form2( Xh, Xh, A ) +=
      on( markedfaces(mesh, mesh->markerName("Bord")) , u1, B, cst(0.) );

 //----------------------------------------------------------------------------------------------------------//

  backend_ptrtype backend(backend_type::build( this->vm() ));
  vector_ptrtype U( backend->newVector( Xh ) );
  backend->solve( _matrix=A, _solution=U, _rhs=B );
  u1 = *U;

}


}//end namespace schwarz_cas2





/**
 * main function: entry point of the program
 */
int
main( int argc, char** argv )
{


  //dim, N, OrderGeo
  schwarz3d_cas2::Localisation<3,2,1> moncode( argc, argv,
                                               schwarz3d_cas2::makeAbout(),
                                               schwarz3d_cas2::makeOptions() );

  //  moncode.run();
  moncode.runMultiDomain();
  //  moncode.runTest();
}



